Formal Foundations of Dressed Time-Dependent Density-Functional Theory for 

Many-Electron Excitations 



Miquel Huix-Rotllant and Mark E. CasidfQ 
Laboratoire de Chimie Theorique, 
Departement de Chimie Molecularie (DCM, 
UMR CNRS/UJF 5250), 
Institut de Chimie Moleculaire de Grenoble (ICMG, FR2607), 
Universite Joseph Fourier (Grenoble I), 
301 rue de la Chimie, BP 53, 
F-38041 Grenoble Cedex 9, France 
(Dated: August 10, 2010) 

In so far as time-dependent density-functional theory (TDDFT) provides an exact formalism for 
calculating molecular absorption spectra, TDDFT should be able to describe not only 1-electron 
excited states but also 2-electron, 3-electron, etc. excited states which show up in molecular spectra 
by borrowing intensity from 1-electron excited states. This requires going beyond conventional adi- 
abatic TDDFT where the exchange-correlation kernel is frequency independent, /„, to include an 
appropriate frequency dependence, f xc (u). Maitra, Zhang, Cave, and Burke gave a heuristic deriva- 
tion of a frequency-dependent correction to adiabatic TDDFT (an approach they called "dressed 
TDDFT") designed to bring in one double excitation [J. Chem. Phys. 120, 5932 (2004)] and Casida 
showed how dressed TDDFT might be generalized in the form of a polarization propagator (PP) 
correction to adiabatic TDDFT [J. Chem. Phys. 122, 054111 (2005)]. This paper presents a further 
exploration of the PP approach to dressed TDDFT using an approach which is significantly more 
rigorous than previous dressed TDDFT work and previous work on a PP correction to adiabatic 
TDDFT. A link is also made with work based directly upon the Bethe-Salpeter equation. At first 
order we recover the exact exchange result of Gorling [Int. J. Quant. Chem. 69, 265 (1998).] A 
second-order treatment brings in 2-electron excitations. An important result of Gonze and SchefHer 
[Phys. Rev. Lett. 82, 4416 (1999)] emerges as a trivial consequence of this formalism. Concrete 
formulae for the nonadiabatic correction are derived at the level of the second-order polarization 
propagator (SOPPA) and algebraic diagrammatic construction (ADC) approaches. The example of 
butadiene is used to make a connection with the pioneering dressed theory of Maitra et al. 



I. INTRODUCTION 

The time-dependent extension of density-functional 
theory (DFT) is currently providing many opportunities 
and challenges for modeling electronic time-dependent 
properties and excitation spectra. [H-Q The conventional 
DFT of Hohenberg, Kohn, and Sham, though a for- 
mally exact theory, is limited to static ground state prop- 
erties. @, U Time-dependent DFT (TDDFT) is a for- 
mally exact extension of DFT designed to treat time- 
dependent phenomena, @ from which excited states may 
be extracted using linear response theory. Q Linear- 
response TDDFT (LR-TDDFT) calculations have be- 
come a method of choice for calculating the spectra 
of medium and large sized molecules. Limitations on 
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TDDFT arise because the exact exchange-correlation 
(xc) functional must be approximated in practice. Thus 
conventional LR-TDDFT works best for low-energy 1- 
electron excited states which do not involve too much 
charge transfer and are not too delocalized in space. [||,|9| 
This paper is a contribution towards overcoming the 
restriction to 1-electron excitations. In particular we 
propose a method for placing the "dressed TDDFT" 
of Maitra, Zhang, Cave, and Burke within the 
well-defined context of many-body perturbation theory 
(MBPT). After some well-defined approximations, we 
will derive an expression for the frequency-dependent xc- 
kernel, f xc (u), of LR-TDDFT at the level of the second- 
order polarization propagator (SOPPA) jllT4l4l and al- 
gebraic diagrammatic construction (ADC) [15], |T(| ap- 
proaches. The result is a powerful formalism which en- 
globes several previous results from the literature but 
which also allows us to go beyond these results to arrive 
at a rigorous and quite general polarization propagator 
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formulation of dressed TDDFT. 

In the most applications of TDDFT, we first begin with 
a molecule in its ground state where it is described by the 
usual Kohn-Sham equation, Q 



h H + Ware (1) J 1pr(X) = £rVv(l) , 



(1.1) 



and introduce an external time-dependent perturbation, 
Uappi(l)- The time-dependent Kohn-Sham equation Q is 
then. 



hH + V xc (l) + V app l(l)) 4> r (l) = i 



.9^(1) 



dt 



(1.2) 



[We use Hartree atomic units, fi = m e = e = 1, and 
space-spin coordinates throughout. Also i = (r^Ci) rep- 
resents the space and spin coordinates of particle i while 
i = (i,U) (bold face) is the corresponding space-time 
coordinate.] Practical applications of TDDFT almost al- 
ways make use of the adiabatic approximation (AA), 



,AA 



(1) 



5E xc [p t ] 



(1.3) 



and pt(x) is the density p(x, t) considered as a function 
only of x. This amounts to assuming that the response 
of the xc-potential, v xc , to any temporal change in the 
density, is instantaneous and without memory. The usual 
form of LR-TDDFT Q involves solving a pseudoeigen- 
value problem, 



A(w) B(w) 
B*(w) A*(w) 



1 
-1 



(1.4) 



where, 



A ia j b (u}) = <i.. />'„.<.<„., + (ia\f H xc(u)\jb) 
B ia , b:j {uj) = (ia\f Hxc (uj)\bj) . (1.5) 

Here, }hxc{u) = fn + fxc(u) where /n(l,2) = l/r 12 is 
the Hartree kernel and f xc {w) is the Fourier transform of 
the xc-kernel, 



/ a , c (l,2;ti-t 2 ) 



5p{2) 



(1.6) 



Also we have introduced the compact notation, 

e rs ... , uv ... = (e r + e s H ) - (e u + e v H ) . (1.7) 

Equation (|1.4[) has paired excitation and de-excitation 
solutions. Its eigenvalues are (de-)excitation energies (or 
frequencies since Ti = 1), 



UJI 



Ei-Eq 



(1.8) 



and the vectors X and Y provide information about tran- 
sition moments. In particular the oscillator strength, 



/^-^KOlrl/)! 2 , 



(1.9) 



of the transition with excitation energy wj may be cal- 
culated from X/ and Y/. p} When the AA is made, the 
A and B matrices become independent of frequency, and 
it is easy to show that the number of solutions is equal 
to the number of 1-electron excitations, albeit dressed to 
include electron correlation effects. Allowing the A and 
B matrices to have a frequency dependence allows the 
explicit inclusion of 2-electron excited states. 

This is easiest to understand within the Tamm-Dancoff 
approximation (TDA). The usual AA TDA equation, 



AX = wX, 



(1.10) 



resembles a configuration interaction (CI) equation, [Vi 



(H-£ l)C = wC 



(1.11) 



However, while Eq. (|1.10[) is restricted to single excita- 
tions, Eq. (jl.lip includes not only single but also higher 
excitations. Nevertheless there is a simple transforma- 
tion which recasts Eq. (jl.lip into the form of Eq. (jl.101) , 
but with a frequency-dependent A(w). Thus the matrix 
A CI = H — EqI may be partitioned into 1-electron ex- 
citations ("singles", denoted here by the subscript "1") 
and other excitations (typically "doubles", but possibly 
also higher excited terms, denoted here by "2+") as, 



A CI 



A ci 

\ CI 

A 2+,2+ 



C2+ 



Ci 

C2+ 



(1.12) 



provided we can ignore any coupling between the ground 
state and our singles, doubles, and higher excited states. 
According to standard Lowdin (Feshbach) partitioning 
theory (see, for example, Ref. fl8j). 



* CI , A CI 



1 ci 



1 A CI 



(wl 2 +,2+ - A 2+j2 

It is tempting to identify, [13] 

A(w) w Af [ + Af i + (usl 2+ 2+ - 2 

(1.14) 

This is effectively what Maitra, Zhang, Cave, and 
Burke did in their dressed TDDFT, except that they 
went further and considered only the case of coupling 
between a single 1-electron excitation (S) and a single 



Ci =uCi. 
(1.13) 



kCI 

A 2+,l 



3 



2-electron excitation (D). This is the model shown 
in Fig. HI Let us try to understand the basic physics. A 
bright singly-excited state with excitation energy wg and 
oscillator strength fs = 1 interacts with a dark doubly- 
excited state with excitation energy lod and oscillator 
strength fjj = via a coupling matrix element x. The 
CI problem is simply, 



uj s x 
x ujo 



C s 
C D 



C s 
C D 



UJ S 


X 


( cos 9 




f cos 9 \ 






) = w.| 




X 


UJ D 


\ sind 


v sin 9 ) 



It has two solutions, namely, 



and 



Backsolving we find that, 



UJ S 


X 


( -sine ' 




-sin 9 \ 










X 


UJ D 


^ cos 9 


cos 9 J 



lus = ui a cos 9 + ujb sin 



u>e> = uj a sin 



uib cos 



(1.15) 



(1.16) 



(1.17) 



(1.18) 



which establishes the unperturbed quantities in terms of 
the solutions of the perturbed problem. We see the av- 
erage excitation energy is conserved in the coupled prob- 
lem, 



uj a + uj b = uj s + u D 



(1.19) 



Something similar occurs with the oscillator strengths. It 
is easy to show that, 



/ a = - Wa |(0|r|a)| 2 



^ a |(O|r|S)| 2 cos 2 



fb 

Summing gives 



|r|S)| 2 sin 2 0. (1.20) 



fa + fb = g K COS 2 I 



u)b snr 



9) |(0|r|5)| s 
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W5)| 2 =/ s . 



(1.21) 



Hence oscillator strength is also conserved. This leads to 
the common interpretation that the coupling "shatters 
the singly-excited peak into two satellite peaks." 

Now let us try to obtain the dressed TDDFT form for 
the nonadiabatic correction to the exchange-correlation 
kernel. Equation (|1 . 13|) reads, 



Cs — ujCs ■ 



(1.22) 



and u) a and w;, are thus solutions of the equation, 

x 2 

LO = bJ S H . 

UJ — Ul£, 



(1.23) 



Comparing with the diagonal TDA LR-TDDFT within 
the two-orbital model, 



w = e<M + (ia\f H xc(u)\ia) 



shows that, 



(ia\f Hxc (uj)\ia) = (uj s - e 0)i ) + 



UJ — U)Q 



(1.24) 



(1.25) 



Maitra et al. interpreted the first term as the adiabatic 
part, 



JHxc — - ta.i , 

and second term as the nonadiabatic correction, 

,.2 



Since 

x = (uj a 
it is easy to show that, 



UJ — UJQ 



UJb) cos ( 



UJ S UJ D - UJ a UJ b . 



(1.26) 



(1.27) 



(1.28) 



(1.29) 



which is the form of the numerator used by Maitra et 
al. Equation (|1.27[) then provides a model for including 
a single double excitation into conventional TDDFT cal- 
culations. This basic idea has now been used by several 

groups. 0, mm 

Maitra et al. recognized that dressed TDDFT should 
ultimately be formulated in terms of linear response the- 
ory rather than the wave function CI method, so they 
gave their justification in terms of generalized suscepti- 
bilities. However their formula is the same as the one pre- 
sented above and apparently involves a double-excitation 
energy, <jd, obtained from wave function theory. Casida 
recognized that a correct formulation of a nonadiabatic 
correction term to adiabatic LR-TDDFT should come 
from linear response theory and so used the polarization 
propagator formalism to devise a polarization propagator 
correction. [22| This earlier theory had some difficulties 
which arc overcome in the present formulation. 

The approach taken here has been termed " ab initio 
DFT" by Bartlett. @, 0] At the exchange-only level 
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FIG. 1. Two- level model used by Maitra et al. in their ad hoc 
derivation of dressed TDDFT. See explanation in text. 
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the terms optimized effective potential (OEP) |25|, [26( 
or exact exchange (EXX) [5J [28| are also used, and 
OEP is sometimes also used to include the correlated 
case. [2!| [3(| The objective is to derive expressions for 
DFT quantities, such as the xc-potential, v xc , and the 
xc-kernel, f xc [uj), from many-body perturbation theory 
(MBPT). (In contrast to the way the term is used in 
solid-state physics, Bartlett uses the term "ab initio" as 
it is normally used in quantum chemistry to make a dis- 
tinction between DFT and typically Hartree-Fock-based 
first-principles theories.) This approach does not give 
functionals of the density, but can give useful insight into, 
and computational checks of, the behavior of DFT quan- 
tities. We will assume that the the xc-potential is known 
and concentrate on finding MBPT expressions for the 
frequency-dependent xc-kernel. 

Previous work along these lines has been carried out for 
the static kernel by directly taking the derivative of the 
OEP energy expression with the constraint that the or- 
bitals come from a local potential. This was first done by 
Gorling in 1998 |3l| for the full time-dependent exchange- 
only problem. In 2002, Hirata et al. redid the derivation 
for the static case. [32| Later, in 2006, a diagrammatic 
derivation of the static result was given by Bokhan and 
Bartlett, [24[ and the functional derivative of the kernel, 
g x , has been treated by Bokhan and Bartlett in the static 
exchange-only case. [33| 

Our approach is more closely related to a second way 
to derive the same equations which may have the advan- 
tage of being more transparent and in any event of hav- 
ing a more direct relation with recent developments in 



solid-state theory. [34H39I ] The essential idea starts from 
the fact that f X c(u) satisfies a 2-point Bethe-Salpeter-like 
equation which must be the result of localizing the corre- 
sponding "proper particle-hole self-energy" (also known 
as the "proper particle-hole scattering kernel" and the 
"compact vertex part" [66]) that appears in the 4-point 
Bethe-Salpeter equation in ordinary many-body theory. 
Tokatly, Stubner, and Pankratov have used this concept 
to develop a diagrammatic expansion of f xc (ijj). (4(1 |4lj 
This approach has been extended to include 2-electron 
excitations in recent work by Romaniello et al. using a 
screened interaction approach inherited from solid-state 
polarization propagators. (42j Their approach is not en- 
tirely satisfactory as they show that it gives rise to extra 
unphysical solutions, a problem which we believe may 
arise from an imbalanced treatment of antisymmetry due 
to an asymmetric treatment of direct and exchange dia- 
grams when using a screened interaction. This belief is 
confirmed by recent work of Sangalli, Romaniello, Colo, 
Marini, and Onida. (43|. 

The basic idea of our approach, which does balance 
contributions from direct and exchange diagrams, is de- 
scribed in the next section. We do not in fact local- 
ize the Bethe-Salpeter equation (BSE) though we do 
the equivalent. The problem is that we wish to work 
in the frequency representation and that the true BSE 
(before approximations) is an equation in the time rep- 
resentation whose Fourier transform into the frequency 
representation is fraught with problems. Instead we fol- 
low the superoperator polarization propagator approach 
(SPPA) of Nielsen, J0rgensen, and Oddershede llMl 
and the algebraic diagrammatic construction (ADC) ap- 
proach of Schirmer [la . [l6| and begin with a formally ex- 
act form of the polarization propagator (Sec. IIVH This 
gives a matrix equation whose elements are then deduced 
in an order-consistent manner by expanding and match- 
ing against Feynman diagrams (Se c . IIV1 ) (See also ear- 
lier work by Paldus and Cizek. [44|, |45[) In effect we are 
resuming these Feynman diagrams in a particular form. 
The result is a rcnormalizcd theory which has proven 
itself for quantum chemical applications but which so 
far seems to defy diagrammatic representation. Unfor- 
tunately the second-order SPPA and ADC approaches 
(termed SOPPA and ADC(2) respectively) do not lend 
themselves to solid-state applications where it is likely 
that screening must be explicitly taken into account. 

Our contribution is to generalize and to adapt these 
equations for use in ab initio DFT. Since we do make 
some well-defined (though not strictly necessary) approx- 
imations to keep future computations tractable, it is im- 
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portant to show that these equations are at least as good 
as those previously presented in the literature. This we 
do by indicating how previous literature results emerge 
in a natural fashion. However our equations are very 
general and go beyond previous work. 

Finally we make the connection with dressed TDDFT 
(Sec. EJ. The result is what we believe to be a coher- 
ent generalization of dressed TDDFT to include an in 
principle infinite number of 1-clcctron, 2-electron, and 
higher-electron excited states and which we are tempted 
to call "fully-dressed TDDFT" because it includes sig- 
nificantly more states than the original dressed TDDFT. 
To make the relation between dressed TDDFT and the 
present theory more clear we present the example of bu- 
tadiene, previously treated by Maitra et al, thus facili- 
tating comparison with their pioneering work (Sec. fVf . 
Section I VII contains our concluding discussion. Appen- 
dices have been included where deemed necessary to keep 
the present work self contained (Appendices [Al and |C|) or 
where additional details seem useful (Appendix [Bj. 



II. SPACE, TIME, AND THE EXACT KERNEL 

This section elaborates on the polarization propaga- 
tor (PP) approach to calculating the xc- kernel, initially 
suggested by one of us. Q Not only do we wish to clar- 
ify certain approximations implicit in this previous work, 
but we wish to bring into sharp relief the relation with 
the Bethe-Salpeter equation (BSE) approach to calculat- 
ing the xc-kernel popular in the solid-state physics com- 
munity. [34M39I . |42| | Thus both approaches involve spatial 
localization of nonlocal operators, some of which (vide in- 
fra) results in new frequency dependencies. However, al- 
though the two many-body perturbation theory (MBPT) 
approaches to ab initio DFT are formally equivalent, dif- 
ferences emerge because the BSE approach emphasizes 
the time representation while the PP approach empha- 
sizes the frequency representation. This can and typically 
does lead to different approximations. The discussion of 
dressed TDDFT in the previous section suggests that it 
will be easier to derive pole structure-conserving approx- 
imations needed for treating 2-electron and higher exci- 
tations in the frequency representation than in the time 
representation. This and prior experience with the PP 
approach in the quantum chemistry community (liTfl6j 
have lead us to favor the PP approach. 

The common denominator to all approaches is the fun- 
damental equation of LR-TDDFT which simply says that 
the response of the density of the interacting system to an 



applied perturbation must be the same as the response of 
the density of the Kolm-Sham fictitious system of non- 
interacting electrons to an effective potential involving 
the response of the Kohn-Sham single-particle potential, 
v s . This response is conveniently described in terms of 
susceptibilities. The susceptibility 



X(l,2) 



Sv app i(2) 



(2.1) 



describes the response of the real electrons to the applied 
perturbation 5v app i, 



S P (1) = J X (l,2)5v appl (2)d2. 



(2.2) 



The response of the density of the Kohn-Sham fictitious 
system of nonintcracting electrons is identical but the 
potential is now the Kohn-Sham single-particle potential, 



5p(l) = J Xs (1,2)^(2)^2. 



(2.3) 



In contrast to the interacting susceptibility of Eq. (|2.ip , 
the noninteracting susceptibility, 



X-(l,2) 



Sv s {2) 



(2.4) 



is known exactly from MBPT. Of course the effective 
potential is the sum of the applied potential and the po- 
tential due to the response of the self-consistent field, 

VHxc, 

Sv s (l) = Sv appl (l) + J f Hxc (l, 2)5p(2) d2 , (2.5) 

where Jhxc{^- 1 2) = 6vu xc (l)/Sp(2) is the functional 
derivative of the Hartree plus exchange-correlation self- 
consistent field. Manipulating these equations is facili- 
tated by a matrix representation in which the integration 
is interpreted as a sum over a continuous index. Thus, 

Sp = X5v app l = Xs {Sv app i + fHxcSp) , (2.6) 

is easily manipulated to give a Bethc-Salpctcr-likc equa- 
tion, 



X 7Cs ~\~ XsSiixcX j 

or, written out more explicitly, 



(2.7) 



X (l,4)=x.(l,4) + J x.(l J 2)/ Hxc (2 J 3)x(3,4)d2d3. 



(2.8) 
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We will be making extensive use of MBPT. Perhaps 
the most common and arguably the most basic quantity 
in MBPT is the 1-electron Green's function defined by, 

»G(1,2) = <0|T{V^(1)V^(2)}|0>. (2.9) 

It hardly needs to be defined except to introduce our no- 
tation. Here, the subscript H indicates that the field op- 
erators are understood to be in the Heisenberg represen- 
tation. Also T is the usual time-ordering operator, which 
includes anticommutation in our case (i.e., for fermions), 

rftMl)$r(2)} = 0(h - t 2 )^ H (l)^ H (2) 
-0(*a-ti)$r( 2 )lMl)- 

(2.10) 

The usual MBPT approach to evaluating the suscep- 
tibility, x, uses the fact that it is the retarded form, 

ix(l, 2) = 8{h - i 2 )(0|[pH(l),^(2)]|0) , (2.11) 

of the time-ordered correlation function, 

» X (l,2) = (0[T{pH(l)Pir(2)}|0) > (2.12) 

where, 

p H (l) = ft H (l)$ H (l) - (01^(1)^(1)10) , (2.13) 

is the density fluctuation operator. (See for example 
Ref. H| pp. 172-175 and p. 151.) 

In practice it is much more convenient to work with 
time-ordered, rather than retarded, quantities. For this 
reason, we will generally work with time- ordered Green 
functions, even when our notation may suggest the cor- 
responding retarded quantity. The retarded version may 
be obtained from the time-ordered version of x by simply 
changing the sign of the imaginary part in the frequency 
representation. Since we typically work with finite ba- 
sis sets and real functions, there is no imaginary part in 
our calculations and so the problem does not arise ex- 
cept in principle. However it is best to keep in mind 
the distinction between different types of Green's func- 
tion time-ordering. (See for example the discussion of 
Lehmann representations on pp. 72-79 of Ref. [46[ for 
the distinction between retarded, advanced, and time- 
ordered Green's functions.) 

The BSE-like form of the basic equation ([277]) of LR- 
TDDFT suggests that the most direct link with MBPT 
should be through the BSE itself. Although the MBPT 
literature refers to several related equations as BSEs, the 



original BSE [ijj is probably that of the particle-hole 
(ph) propagator, [lj| which we choose to write as, 

iL(l, 2; 3, 4) = <0|T{7(1, 2) 7 (4, 3)}|0) , (2.14) 

where, 

7 (1, 2) = ^(2)^(1) - <0|T{$r(2)iMl)}|0) , (2.15) 

is a sort of density matrix fluctuation operator (or would 
be if we constrained t\ = ti and t% = t^). 

[It may be useful to try to place L in the context 
of other 2-electron propagators: The 2-electron Green's 
function is (see p. 116 of Ref. |46|). 

G(l, 2; 3, 4) = (-i) 2 (0\T{iP ff (l)i> H (2)^(4)^(3)110). 

(2.16) 

The particle-hole response function, [15| 

i?(l,2;3,4) = G(l,2;3,4) - G(l, 3)G(2, 4) . (2.17) 

Then L is related to R by the relation, 

L(l,2;3,4) = iR(l,4;2,3).] (2.18) 

By construction, 

X (l,2)=i(l,l+;2+,2), (2.19) 

where i + is infinitesimally later than i. In the jargon of 
MBPT, x(l, 2) is a 2-point quantity, while L(l, 2; 3, 4) is 
a 4-point quantity. The space-time analogue of the Har- 
riman contraction operator (Appendix^]), Yst, allows 
us to rewrite Eq. (|2.19[) as a concise matrix equation, 

x = r ST Lrl. T . (2.20) 

Much of the development of the BSE approach to cal- 
culating the xc-kernel done in the solid-state physics 
community (334391 l42l [43| makes use of Hedin's equa- 
tions. |48| These cannot be used directly to determine 
the xc-kernel without at least some manipulation which 
would take us too far afield. Suffice it so say that an 
appropriate BSE for 

L(l,2;7,8) = L S (1,2;7,8) 

+ J L a (l,2;3,4)B Bxc (3,4;5,6)L(5,6;7,8)d3d4d6d6 J 

(2.21) 

or 

L = L s + Ls&HxcL , (2.22) 
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in matrix notation. Here 



iL.(l,2;3,4) = G a (l,4)G,(2,3) 



(2.23) 



is the ph-propagator for the Kohn-Sham fictitious sys- 
tem of noninteracting electrons, and the 4-point quan- 
tity, Sjj xc , may be deduced from a Feynman diagram 
expansion as the proper part of the ph-response function 
"self-energy" . 

Combining Eqs. (|27f]l . (|2~2"0"j) . and (j2~2^|) gives the for- 
mally exact equation, 



TsrLTl 



- )fxc[ i-ST^ 



YsTL s a xc LY 



ST 



(2.24) 



Note that the Hartree part of Sh X c is just Yjy T /#YsT 
and so cancels with the Hartree part of Y^ ST fuxc"^ST- 
Of course, even ab initio DFT involves approximations 
to truncate infinite sums and otherwise keep calculations 
tractable. A common approximation is to replace the 
computationally difficult many-body quantity L with the 
computationally much simpler single-body quantity L s 
to obtain, 



ST 



~£sTL s E xc L s 'Tg T . (2.25) 



This is still very complex because H xc , is a 4-time quan- 
tity (or 3-frequency quantity after Fourier transforma- 
tion [iH). Equation (|2.25p is closely related to the 
"Nanoquanta approximation" (NQA), 

stLsY^j^J f xc (y S tL s Y\, t 

i ~5G~ ) Lsr s T > 

(2.26) 

so named by Lucia Reining because it was simultaneously 
derived by several different people [334361 [H, [H, |49| in- 
volved with the Nanoquanta group. The NQA simplifies 
still further in the usual GW approximation for the self- 
energy £ xc . Although the NQA has proven successful for 
extended systems, it is not at all clear that it is appropri- 
ate for finite systems where vertex effects are likely to be 
important. However its most severe drawback, and the 
one which prevents us from using it here, is that it is ill 
adapted for treating 2-electron and higher excitations. 

We have chosen to seek more appropriate approxima- 
tions within the PP approach. This approach is based 



upon the 2-time quantity, 

n(l, 2; 3,4; t - t') = L(lt, It; it' , At') , 



(2.27) 



which makes it impossible from the beginning to write 
down a BSE for the PP (though Bcthc-Salpctcr-like equa- 
tions do arise within the PP approach.) Written out ex- 
plicitly, 

TI(l,2;3,4;i-t') 
= (0|T{^(2t+)^(lt)^(3t' + )^(4t')}|0) 

- (o|r{^(2t+)^(it)}|o)(o|r{v^(3t /+ )^(4t')}|o) . 

(2.28) 

[The second term is often dropped in the definition of 
the PP. It is there to remove u> — excitations in the 
Lehmann representation. (See for example pp. 559-560 of 
Ref. [46|.)] The retarded version of the PP is the suscep- 
tibility describing the response of the 1-electron density 
matrix, 

7(l,2;t) = (0$t(2f$(lt)|0), (2-29) 
to a general (not necessarily local) applied perturbation, 

*7(l,2;t) 



n(l,2;3,4;<-t') - 

dw app i(?>,A\t') 

which is a convolution. After Fourier transform, 



(2.30) 



<5 7 (l,2;u;)= / n(l, 2; 3, 4; u))5w app i(3, 4; u) d3d4 , 

(2.31) 



or. 



6f{w) = U(u)Sw app i(uj) , 



(2.32) 



in matrix form. The connection with with DFT is made 
by recognizing that the diagonal elements of the density 
matrix are the density, 



p(l;w) = 7 (l ) l;w) 
p{u) = T7(w) , 



(2.33) 



and DFT is typically only concerned with local operators, 

Sw app i(l, = 6(1 - l')Sv app i(l;uj) 

5w ap pi(u) = Y^Sv app i(uj) . (2.34) 

Here we have introduced the original form of the very 
convenient collapse and expansion operators of Harriman 
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(Appendix [A} . Applying this information to Eq. (|2.31[) 
gives that, 



so that, 



xM = Tn(w)rt 



(2.35) 



The basic LR-TDDFT Eq. {27]) is a convolution which 
when Fourier transformed gives, 

X(l,4;w) = x»(l,4;w) 

+ /" x a (l,2;c;)/ff 1BC (2 > 3;«)x(3 J 4;w)d2d3 J 

(2.36) 

or in matrix form, 

xM = x„(w) + Xs(u)fHxc[u)x(u) ■ (2-37) 



Applying Eqs. (|2.33[) and (|2.34[) leads to the exact equa- 
tion, 

(Tn,(w)Tt)/ Hac ( w )(Tn(w)Tt) 

= T(n(w) -II,(w))T t .(2.38) 

Further approximations are useful in order to minimize 
the number of times the computationally difficult full PP 
appears and to facilitate interpretation. However the pre- 
cise nature of these approximations typically differ de- 
pending upon whether the interesting frequencies are far 
from any higher-order excitations (off-resonant regime) 
or near a higher-order excitation (resonant regime). In 
discussing these two cases, it is helpful to rewrite the 
exact case as, 



Exact: 

f Hxc (uj) = a s (uj) (njV) - n-V)) At( w ) , 



where the localizer, 

A(w) = (Tn(w)Tt) -1 Thft 



(2.39) 



(2.40) 



Note that this only makes sense if we take care to elim- 
inate the null spaces of Yn s (ci;)Yt and of YII(w)Yt 
before taking the inverse (i.e., singular value decomposi- 
tion.) See Sec. lIIII and Appendices lAl and [51 for a detailed 
discussion of the localizer. 

Optimized effective potential theories typically make 
use of the Born approximation, 

xM = XsM + Xs{u)f H xc{u)Xs{u) , (2.41) 



/hic(w) = (3C»(w)) 1 (xH - XsM) (XsH) 1 



This is our first approximation, 



(2.42) 



1st Approx.: 

f Hxc (oj) = a s (lu) (nrV) - n- 1 ^)) a\ /2 (lu) 
= (rn a (w)Tt) -1 (n( w ) - n a ( w )) (Yn s ( w )Yt)" 1 . 

(2.43) 

where only one II has been replaced with II S in the lo- 
calizer, 



A 1/2 (u) = (Tn < (w)T t ) TII(w). 



(2.44) 



It is expected to work well in the off-resonant regime. 
As we shall see, this first approximation gives Gorling's 
exact exchange (EXX) kernel for TDDFT. HH On the 
other hand, the poles of the kernel in this approxima- 
tion are a priori the poles of the exact and independent 
particle PPs — that is, the true and single-particle excita- 
tion energies — unless well-balanced approximations lead 
to fortuitous cancellations. The discussion of dressed 
TDDFT in the introduction shows that this is unlikely 
to lead to accurate 2-electron and higher excitations. 

To correct the poles we make an additional structure- 
conserving approximation, equivalent to the PP Born ap- 
proximation, 

n(w) = II a (w) + U s {oj)K Hxc {cu)U s (lj) 

K Hxc = UJ 1 {u)-U- 1 (u). (2.45) 

This is our second approximation, 



2nd Approx.: 



n-V)) At(o 



(2.46) 



Not only does the static localizer enter into the static 
OEP problem (Appendix \K§ , but the dynamic localizer 
also enters naturally into the time-dependent OEP prob- 
lem when the initial state is the static ground state. [3l[ 
Equation (|2.46[) simply reads that /ij ac (a;) is a spatially 
localized form of the quantity, Kh xc (ui) — nj 1 (o;) — 
II~ 1 (a;). This is nothing but the PP analogue of the ba- 
sic approximation (|2.25[) used in the BSE approach on 
the way to the Nanoquanta approximation. 
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III. LOCALIZATION 

At the risk of repetition, let us recapitulate what we 
have done: We have made two Born approximations, 

n(w) = n.(w) + u s {oj)K Hxc {cu)n s (Lj) 

XM = X*(u) + Xs{u)f H xc{u)Xs{u) , (3.1) 
and have thereby arrived at the equation, 

/jtscM = A s (uj)K Hx c(u)Al(u>) . (3.2) 

The significance of this equation is that the ab ini- 
tio many-body problem has now been separated into a 
purely PP part represented by ~K.Hxc(^>) and a frequency- 
dependent localization part represented by A s (uj) which 
allows us to recover a TDDFT quantity from ab ini- 
tio theory. This makes manifest earlier ideas suggest- 
ing that localizing in space is necessarily associated with 
"derealization in time" (i.e., additional frequency depen- 
dence. [2!|[50| Since H(u) and x( w ) have the same poles 
and the same transition densities (but not the same tran- 
sition density matrices!), spatial localization is not neces- 
sary if our only goal is to obtain excitation energies and 
transition densities beyond the TDDFT AA, since it then 
suffices to carry out an ordinary DFT-based PP calcula- 
tion. Rather the extra localization step is needed if we 
wish to remain within TDDFT and hence study intrinsic 
TDDFT quantities such as the xc-kernel. In such a case, 
exploration of the behavior of the dynamic localizer is 
essential. In this section, we give the explicit form of the 
dynamic localizer from which we easily deduce an impor- 
tant relation first noted by Gonze and Scheffler [5l| . This 
Gonze-Scheffler relation was used in the original justifi- 
cation of dressed TDDFT by Maitra, Zhang, Cave, and 
Burke. Its importance to the formal foundation of 
dressed TDDFT is re-examined in Sec. |V] Eventually it 
will be desirable to go beyond dressed TDDFT to carry 
out finite-basis calculations with the dynamic localizer 
(not done in this paper) and this and a static approxi- 
mation are discussed in Appendix |B| 

The static localizer is intimately related to the static 
OEP approach and has been discussed in Appendix [3] 
The dynamic localizer, 



A(w) = (Tn(w)T t ) TU(u) 



(3.3) 



arises quite naturally in the context of the time- 
dependent OEP problem. According to Runge-Gross 
(RG) theory, Q the exact time-dependent xc-potential, 
v X c{t), is not only a functional of the density, p(t), but 



also of an initial condition which can be taken as the 
wave function, ^(to), at some prior time to. On the 
other hand, linear response theory begins with the static 
ground state case where the first Hohcnbcrg-Kohn theo- 
rem tells us that the wave function is a functional of the 
density ^(io) = ^[pt }- Gorling has pointed out that this 
greatly simplifies the problem 31| because we can then 
show that, 

= yn.(l,l;2,3;w)E x (2,3)d2d3. (3.4) 

We can write this equation in matrix form using the Har- 
riman's collapse operator as 



Tn s ( w )rtv I (w) = rn 8 ( w )s x , 



(3.5) 



Clearly the dynamic exchange potential is the Hartree- 
Fock exchange operator localized by the dynamic local- 
ization operator, 



v x (w) = A s (o;)S a 



(3.6) 



For the case of the non-interacting susceptibility, we 
can easily derive an expression for the dynamic localizer. 
Since 



n s (l,2;3,4; W ) = ^^ 

i a 
occ virt 

-EE 



^(1)^(2)^(3)^(4) 



^(1)^(2)^(3)^(4) 



(3.7) 



we can express the kernel of YII s (a;) as 

occ virt 



(TII.)(l;2,3;a;)=X;£ 



^(1)^(1)^(2)^(3) 



^(1)^(1)^(2)^(3) 



w + e a i 



(3.8) 



Also, the kernel of TII s (w)T^ is just, 



(TILTt) (1;2;oj)=Y^J2 

i a 
occ virt 

-EE 



^(1)^(1)^(2)^(2) 



^(1)^(1)^(2)^(2) 



(3-9) 
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Like the susceptibility, The two operators have poles at 
the independent particle excitation energies to = ±e a ,i = 
±(e a ~ e-t)- 

(Here and throughout, unless otherwise stated, we use 
the orbital index convention, 

a, b, c, • • • ,hi,j,k,l,m,no,p,q,--' ,z, (3.10) 

unoccupied occupied free 

where free orbitals may be either occupied or unoccu- 
pied.) 

In order to construct the dynamic localizer, the ker- 
nel (|3.9p has to be inverted. This is not generally possible 
to do this analytically, though it can be done in a finite- 
basis representation (Appendix [Bj . However Gonze and 
Schcfiler have noted that exact inversion is possible in 
the special case of a frequency, to = e b _j, of a pole well 
separated from the other poles. [5l[ Near this pole, the 
kernels, TII s (w) and TIl s (w)T^ are each dominated by 
single terms, 

^(1)^(1)^(2)^(3) 

(Tn s ) « J - 

U> - €b,j 

(Tn.Tt) (1; 2; u) » 3 b 3 . (3.11) 

CO - €b,j 

Thus Eq. (|3.5[) becomes, 

!MMW WsK(e(j)WjH MM£) W6 |f ; , Wi) , 

tO — €b,j tO — Eb.j 

(3.12) 

with the approximation becoming increasingly exact as 
to approaches eb t j. Hence, 

(^b\v x (eb,j)\ipj) = (fpbl^xlipj) • (3.13) 

More generally for an arbitrary dynamic kernel, 
K(l,2;u>), 

(^*|A(e bj )X(e hj )) = tyjiKMtyt) , (3.14) 

and we can do the same for — €b t j, obtaining 

(^r b \K~^i)K{-e b . 3 )) = (^\K(-e M )\A) (3.15) 

We refer to these last two equations as Gonzc-Schcfflcr 
(GS) relations, since they were first derived by these au- 
thors. [Hl| These relations show that the dynamic local- 
izer, A s (lo), is pole free if the excitation energies, e a< i, 
are discrete and nondegencrate and suggests that the 
dynamic localizer maybe a smoother function of to than 



might at first be suspected. Equation (|3.13[) also shows 
the importance of the dynamic character of the localizer 
as this is what allows a local operator at certain special 
frequencies to have exactly the same matrix elements as 
the original nonlocal operator which was being localized. 

We can now return to a particular aspect of Casida's 
original PP approach [22| which was failure to take 
proper account of the localizer. The importance of the 
localizer is made particularly clear by the GS relations in 
the case of charge transfer excitations. The single-pole 
approximation to the i — > a excitation energy is, 

to = e a ,i + (ia\A(e ai i)K xc (e a j)A' i (e a j)\ai) 

= e<M + (aain; 1 ^) - n- 1 ^)!**) . (3.16) 

Had the localizer been neglected, then we would have 
found incorrectly that, 

w = e-a.i + {ia\U~ 1 (e m ) - n _1 (e a!l )|ai) . (3.17) 

While the latter reduces to just e a i for charge transfer 
excitations at a distance (because ipiijj a — 0), the former 
does not. [52[ However, for most excitations the overlap is 
non-zero. In such cases and around a pole well-separated 
pole the localizer can be completely neglected. 

Away from the poles, the behavior of the localizer be- 
comes much more complicated. It is evident that the 
range space of the operators TII s (w)T' and ~TH s (to) 
and of the localizer A s (lo) is the product basis ipiipa made 
by multiplying together all occupied, ipi, and all unoccu- 
pied, tp a , orbitals. It can be shown that this product 
basis is complete except for constant functions (because 
/ tp a (l)ipi(l) dl — 0) as long as the orbital basis set is 
complete (Appendix of Ref. [HI). It is, in fact, over- 
complete. Even for finite basis sets, the product basis 
set typically contains linear dependencies. This makes it 
difficult to solve the localization equation 

v(u) = A(w)E(w) , (3.18) 

in full generality because we are really solving a least 
squares problem that seeks to maximize the component 
of v(co) in the relevant range space. (Note that we have 
dropped the subscript x and allowed S to have a fre- 
quency dependence, so that Eq. (|3.18[) is intended to be 
general.) This problem is discussed in more detail in Ap- 
pendix [B] 
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IV. PERTURB ATIVE APPROXIMATIONS TO 
THE POLARIZATION PROPAGATOR 



A. Superoperator Polarization Propagator 
Approach 



In this section we review the MBPT techniques needed 
for the calculation of the key quantities for constructing 
of xc-kernels within the first and second approximations. 
In the case of the first approximation, this means deriving 
an expression for the quantity II (cj) — T1 s (lo). For the sec- 
ond approximation, this means deriving an expression for 
the quantity II^^u;) — n _1 (u;). Although it is clear that 
the zero-order approximation should be the Kohn-Sham 
fictitious system of non-interacting electrons, we have 
still had to make some choices. In particular, we have 
chosen to work with MBPT expansions in terms of the 
bare electron repulsion [or more exactly the "fluctuation 
potential" (vide infra)], rather than the screened interac- 
tion used in solid-state physics. [13, S3] The main advan- 
tage of working with the bare interaction is a balanced 
treatment of direct and exchange diagrams, which we feel 
is likely to be especially important for treating two- and 
higher-electron excitations. While we will automatically 
include what the solid state community refers to as ver- 
tex effects, the disadvantage of our approach is that it 
is likely to break down in solids when screening becomes 
important. The specific approach taken here is the now 
well established second-order polarization propagator ap- 
proximation (SOPPA) of Nielsen, J0rgensen, and Odd- 
ershede jTTMl^ ] . The usual presentation of the SOPPA 
approach is based upon the superoperator equation-of- 



motion (EOM) approach previously used by Casida. [22 1 
However the SOPPA approach is very similar in many 
ways to the second-order algebraic diagrammatic con- 
struction (ADC(2)) approach of Schirmer [ll| [l6| and 
we will not hesitate to refer to this approach as needed 
(particularly with regard to the inclusion of various dia- 
grammatic contributions.) What is new, aside from ex- 
plicit comparisons with appropriate TDDFT literature, 
is the change from a Hartrce-Fock to a Kohn-Sham zero- 
order picture and the concomitant inclusion of (many) 
additional terms. Nevertheless it will be seen that the 
final working expressions are fairly compact. 



The PP is given in a molecular orbital basis as, 
11(1,2, 3, 4; t-t') =J]n, w (i-Oi(lWI(^;(%W 1 

pqrs 

(4.1) 

where 

- n sr>?p (t - = i6(t - t')(0\f\j(t)§H(t)q\j(t')pH(t')\0) 
+ i6(t' - t)(0\qUt')PH(t')rUt)s H (t)\0) . 

(4.2) 

As explained in Ref. [46[ pp. 559-560, this change of 
convention with respect to that of Eq. (|2.28[) turns out 
to be more of a convenience than an inconvenience. Also 
note that, since the PP depends only upon the time dif- 
ference, t — t', we can simplify the above expression with- 
out loss of generality by shifting the origin of the time- 
scale so that t' = 0. We wish to seek an elegant matrix 
formulation for II. 

The first step makes use of superoperators. A super- 
operator X is defined by its action on an (Hubert-space) 
operator A as 



xA = [X,A] =xA-Ax. 



(4.3) 



When X is the Hamiltonian operator, H, one often 
speaks of the Liouvillian. An exception is the identity 
superoperator, 1, whose action is simply given by, 



1A 



(4.4) 



The Hciscnberg form of orbital creation and annihilation 
operators is easily expressed in terms of the Liouvillian 
superoperator, 



PH (t) = e m P e- mt 



e mt p. 



(4.5) 



Then 



TL sr , qp (t) = ie(t)(0\ 

+ i6{-t)(0\tfp 



e iHt (fU) 



e im (ft si 



10} • 



(4.6) 



Taking the Fourier transform (and making use of appro- 
priate convergence factors) gives, 

- n sr ,, P H = (pmiwi + Hy^s), (4.7) 
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where we have introduced the superoperator metric, 

(A\X\B) = (0\[A\[X,B]}\0). (4.8) 
It may be useful to note that, 



n 



sr.qp 



(uj) = n rSjPg (w) 



(4.9) 



follows as an easy consequence of the above definitions. 
Moreover since we typically use real orbitals and a finite 
basis set, the PP is a real symmetric matrix. This allows 
us to simply identify II and the superoperator resolvant, 



n„ 



,(w) = (pTglM + tf)- 1 !^) 



(4.10) 



Since matrix elements of a resolvant superoperator are 
harder to manipulate than resolvants of a superoperator 
matrix, we will transform Eq. (|4.7p into the later form by 
introducing a complete set of excitation operators. The 
complete set 



{Tt} = {T\ ; ; ...} = , V& ; a^j , Vtfb ; ...} 

leads to the resolution of the identity (RI), 
1 = |T t )(T t |T t )- 1 (T t |. 



(4.11) 
(4.12) 



We have defined the operator space differently from the 
previous work of one of us to be more consistent with 
the literature on the field of PP calculations. The differ- 
ence is actually the commutation of two operators which 
introduces one sign change. 

Insertion into Eq. (|4. 7[) and use of the relation, 

(T t |(a;l+F) _1 |T t ) = (T f |T+)(Tt luji+H^)- 1 (T* |T+) 

(4.13) 

then gives, 



n 



sr,qp 



(p f q\T^)(T^\iol + HlT^-^T^f^s) . 

(4.14) 

The resolvent is now defined as a matrix over an under- 
lying configuration space somewhat analogous to what 
is done in the configuration interaction (CI) method fa- 
miliar to quantum chemists [Eq. (|l.lll) ]. However the 
"Cl-expansion" in the EOM method is over an operator 
space, rather than a wave function. 

It is to be emphasized that the PP in Eq. (|4.14[) is exact 
and completely general. Its poles are exact many-electron 
excitation energies and its residues are the corresponding 
oscillator strengths. As such, it must be equivalent to 
the BSE which is often expressed diagrammatically in 
terms of Feynman diagrams. (Diagrams are discussed in 
Appendix O) 



r s 

H 




r s 

IT 


t 




T T 






L 




n 






(T+KW^lTt)- 1 



p q 




FIG. 2. Graphical representation of the localization of the 
4-time propagator of the Bethe-Salpeter equation into 2-time 
quantity of the polarization propagator. See text. 



Figure [2] shows schematically how the PP diagrams are 
related to the BSE diagrams. To pass from the BSE di- 
agrams to the PP diagrams, take the time-unordered ex- 
pansion of the BSE and fix two times so that L becomes 
Id. This is indicated in Fig. [2] by the time-collapse op- 
erator joining the two left most diagrams. Only two 
times are fixed in the second diagram, so it is not yet fully 
time ordered. Taking all possible time orderings joining t 
and t' , including times both before and after these times, 
gives the third diagram in Fig. [2] The result is that some 
parts of the diagrams hang down below the earlier time 
and others are stacked up above the latter time. These 
"dangling boxes" correspond to the terms (p'q\T>) and 
(T'\f's) and represent the initial and final state corre- 
lation effects present in (0| and |0). The multiple time 
ordering also leads to multiple recrossings of the t and t' 
lines. This is why the matrix, 



r- 1 ^) = (T t |wi + i?|T t )- 



(4.15) 



involves many more than four indices and so corresponds 
directly to two-, three-, and higher-electron excitations. 

We have presented this procedure as if we simply ma- 
nipulated Feynman diagrams. In reality we expanded the 
matrices using Wick's theorem with the help of a home- 
made Fortran program. Thee result was a series of 
algebraic expressions which were subsequently analyzed 
by drawing the corresponding Feynman diagrams. This 
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leads to about 200 diagrams which we ultimately resum 
to give a more compact expression. It is the generation 
of this expression that we now wish to discuss. 

Let us analyze this expression for the PP accordin g to 
the order of excitation operator. Following Casida, [22[ 
we partition the space as, 

-IW W ) = (Vg|Tl) (ptfflT^)) T-V) ( ( ^'f f ?i) 

where T 2+ corresponds to the operator space of two- 
electron and higher excitations and 



Tl,2+ 

r 2+ , 2 +M 



has been blocked, 

r w (w) 



(tT| w I + #|t}) 



(4.17) 



(4.18) 



Using the well-known expression for the inverse of a two- 
by-two block matrix allows us to transform Eq. (|4.16|) 
into. 



n. 



0q\T\)~(p1q\T 2+) L 2+a+{ u. 



p- 1 ( w )[(Tl|rts)-r 1)2+ r 
+ (^| T t + )r 2 -_ 1 2+ M(Tt + 

where, 



1 

2+,2 



^2 
M(TL 



PH = r M H 



iY 2 +r 2+2+ 



( W )r 



2+.1 



w)r 3 +,i] 

(4.19) 



(4.20) 



Although Eq. (|4.19p is somewhat complicated, it turns 
out that P(u>) plays much the same role in the smaller 
T\ space that T(u) plays in the full T 1 " 

To see how this comes about, it is necessary to in- 
troduce the concept of order. Our zero-order Hamilto- 
nian is that of the Kohn-Sham fictitious system of non- 
interacting electrons and the perturbation defining the 
order is what is left over. Thus, 



H {0) = h KS 

= V + M XC , 

where V is the fluctuation operator, 



v = -^2 {pq\\ rs )p*r*sQ-'%2 (pr\\rq)pU, 

pqrs pgr 



and. 



(4.21) 



(4.22) 



(4.23) 



Here Ti^ F is the HF exchange operator defined in terms 
of the occupied Kohn-Sham orbitals and the integral, 



(pq\\rs) = J ^(lX(2)i(l-Pi2)^(l)V.(2)dld2, 

(4.24) 

is explicitly antisymmctrized. Notice that order param- 
eters have been put as superscripts in parentheses. 

We can now perform an order-by-order expansion of 
Eq. (|4.19j) . Through second order only the T 2 part of 
T 2+ contributes, so we need not consider higher than 
double excitation operators. However we shall make 
some additional approximations. In particular, we will 
follow the usual practice and drop the last term in 
Eq. (|4.19[) because it contributes only at second order 
(Fig. 10, supporting material [HI) and appears to be 
small when calculating excitation energies and transi- 
tions moments using the Hartree-Fock approximation as 
zero-order. [l5l . [55l458| For response functions like dy- 
namic polarizabilities, their inclusion is more critical, im- 
proving the agreement with experiments. [l2| 

We will also have no need to consider the second term 



(ptg|T+) - (ptg| T t + )r^, 2+ Mr 2+>1 ■ (4-25) 

This is because it only enters into the first approxima- 
tion of the previous section if we go beyond first order 
and is killed by the second approximation of the previous 
section. 

This means that for the purposes of this paper we can 
treat the PP in the present work as given by, 

- U sr , qp {uj) = (ptglT^p-VXTi ■ ( 4 - 26 ) 



Comparing with Eq. (|4.14[) substantiates our earlier claim 
that P(uj) plays the same role in the t\ space that T(uj) 
plays over the full space. 



B. First-Order Exchange-Correlation Kernel 

We now turn to the first-order exchange-correlation 
kernel. Our main motivation here is to verify that we 
obtain the same terms as in exact exchange (EXX) cal- 
culations when we evaluate II — II S . [3ll. l59j Since our 
approach is in some ways more general than previous ap- 
proaches to the EXX kernel, this subsection may also 
provide some new insight into the meaning of the EXX 
equations. 
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Since we are limited to first order, only zero- and first- 
order wave function terms need be considered. This im- 
plies that all the contributions due to the space (the 
space of double- and higher-excitations) are zero and sub- 
stantiates our claim that Eq. (|4.26[) is exact to first-order. 
An ordcr-by-order expansion gives, 

-n££>( w ) = (pH\t{)^p^-\u)(t{\^s)^ 

+ (pt$|Tj)(°)p«'- 1 (a;)(Tt|fts)( ) 

(4.27) 

where, 

-H^» = (ptg| T t)(°)(Tt|c 4 ;l+/i xs |Tt)( ).- 1 (Tt|fts)( ) . 

(4.28) 

The evaluation of each of first-order blocks is straightfor- 
ward using the basic definitions and Wick's theorem. 

Let us first consider the P parts. The zeroth-order 
contribution is, 

P £aH = - ^aWac (4-29) 

^» = 0, (4.30) 
and the first-order contribution gives 

p kc,ia = HIM + M ac 8 ik - M ik 5 ac (4.31) 
Pck,ia = (ci\\ak). (4.32) 

The sum of P' -* + P' 1 -* gives the exact pole structure up 
to first-order in the SOPPA approach. 
The zero-order contribution, 

(ptg|Tl)( ) = (T+|Tt), (4.33) 

and the first-order contributions are given by, 



[G^ItI)]^, = 


M 3 

Oik 


(4.34) 


[(P^lTl)]^, = 


M ic . 
—o k j 


(4.35) 


[(pH\T\)] k % a = 


M ka 

Obc 

£k,a 


(4.36) 




M kbx 

Oca ■ 

£ k ,b 


(4.37) 



The PP II(<j) is now easily constructed by simple 
matrix multiplication according to Eq. (|4.27|) . Apply- 
ing the 1st approximation from Sec. [TT] and expanding 




(a) (b) (c) (d) 





H (n) 

FIG. 3. Diagrammatic representation of n(uj) — II s (w). (a- 
i) involve coupling between the particle-hole space, diagrams 
(g,h,m,n) involve coupling between particle-hole space and 
particle-particle and (i,j,k,l) couples the particle-hole space 
with hole-hole. For the rules of the diagrams, see appendix [Cl 

II s (w) — n(w) through first order allows us to recover 
Gorling's TD-EXX kernel. [3l| The most convenient way 
to do this is to expand pW' -1 using, 

(Tjlwl + HlTjr 1 « (T^l + ^^lTj)- 1 
+ (T\\ul + ff^lT^-^TllFWlT^^lwl + H^lTl)- 1 . 

(4.38) 

The result is represented diagrammatically in Fig. [3j The 
corresponding expressions agree perfectly with the ex- 
panded expressions of the TD-EXX kernel obtained by 
Hirata et al, [59j which are equivalent to the more con- 
densed form given by Gorling. [3l| The diagrammatic 
treatment makes clear the connection with the BSE ap- 
proach. There are in fact just three time-unordered di- 
agrams shown in Fig. 0] whose various time ordcrings 
generate the diagrams in Fig. [3] However the "hanging 
boxes" above and below the central T region now have 
the physical interpretation of initial and final state wave 
function correlation (Fig. [2]). Had we applied the 2nd ap- 
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M 



M 



FIG. 4. Topologically different time-unordered diagrams for 
the First-order contribution to the 4-space 2-time Bethe- 
Salpeter equation for the case of a general potential. In the 
case that M = Vhf the diagrams with M circles should be 
dropped. 



proximation of the previous section, then only diagrams 
(a-f) of Fig. [3]would have survived. 

Use of the Gonzc-Schcfflcr relation (see further Sec. 
then leads to, 



KS\ 



KS 



(a\M xc \a) - (i\M xc \i) + {ai\\ia) 



= taf + (ai\\ia) 



(4.39) 



which is exactly the configuration interaction singles 
(CIS, i.e., TDHF Tamm-Dancoff approximation) expres- 
sion evaluated using Kohn-Sham orbitals. This agrees 
with a previous exact result obtained using Gorling-Levy 
perturbation theory. [U [6l| 



C. Second-Order Exchange-Correlation Kernel 

Having verified some known results, let us go on to do 
the MBPT necessary to obtain the pole structure of the 
xc-kernel through second order in the 2nd approximation. 
That is, we need to evaluate II J 1 (u>) — n _1 (tj) through 
second order in such a way that its pole structure is ev- 
ident. The SOPPA/ADC strategy for this is to make a 
diagrammatic T1 s (lo) — Tl(u) expansion of this quantity 
and then resum the expansion in an order-consistent way 
to have the form 



[H( w )-n(4w + "^ = 

n k k — i 
fe=0 i=0 j=0 

(4.40) 

when the Born approximation is applied to the P(w), in 
the same fashion as in the previous section. The number 
of diagrams contributing to this expansion is large and, 



for the sake of simplicity, we will only give the resumed 
expressions for each block. For a cocmpilation of relevant 
diagrams the reader is refered to the supplementary ma- 
terial. [H3] Evidently, after the calculation of each block 
there will be an additional step matrix inversion in order 
to apply the 2nd approximation to the xc-kernel. 

It should be emphasized that although the treatment 
below may seem simple, application of Wick's theorem is 
complicated and has been carried out using an in-house 
Fortran program written specifically for the purpose. 
The result before resummation is roughly 200 diagrams 
which have been included as supplementary material. 

It can be shown that the operator space may be trun- 
cated without loss of generality in a second-order treat- 
ment to only 1- and 2-electron excitation operators. [HI 
The wavefunction may also be truncated at second-order. 
This truncation breaks the orthonormality of the T| 
space: 



(T\\T\) » (Tt|T{)(°) + (Tj|Tj)« + 



1 

-1 



_ (4.41) 

This complication is dealt with by orthonormalizing our 
operator space. The new operator set expressed in terms 
of the original set contains only second-order corrections, 



[o'!l< 2 > 



hid 



(kd\\lb)(dk\\al) 

^kl,bd e kl,da 



y-v / 1 (md\\jc)(ci\\dm) 



mcd 



^mj,cd^im.cd 



E 



M jd M di \ . .. . 

Cj.dCi.d I 



(4.42) 



(Note that we have used the linked-cluster theorem to 
eliminate contributions from disconnected diagrams. For 
a proof for the EOM of the one- and two-particle the 
Green's function see Ref. [62| ) 

We may now proceed to calculate, 



ni 2 > » = (ptglTt^pW'-VXTtlftj) 



Co) 



+ (ptg|T t 1 )(°)p( 1 )- 1 H(T t 1 |fts)( 1 ) 
+ (ptg|Tl)( 1 )p(°)- 1 H(T t 1 |fts)( 1 ) 

(4.43) 

The only new contributions that arise at this level are 
due to the block P' 2 ', which is given by 



>(2) 



r (2) r (l) r (0) -1, x F (l) 
1 1,1 ~ 1 1,2 A 2,2 \ UJ )L 21 . 



(4.44) 
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(We are anticipating the w-dependence of the various T- 

(2) 

blocks which will be derived below.) Since the block T\ [ 
is affected by the orthonormalization procedure, it may 
be useful to provide a few more details. Expanding order- 
by-order, 





[tI[u\ + h (0 \t\]]\q^) 


+ (0<°> 


[tI[(ji + h(°\t{}}\oW) 


+ (0« 


[t\,[ui + h^,t{}]\o^) 


+ (0<°> 


[T\W,[ujl + m o \T\]}\oM), 


+ (o(°) 


[Tl,[ul + H^,Tn\0^), 


+ (o« 


[tUhW,t\]\oW) 


+ (0<°> 


[Tl[HW,T\]\0M) , 



(4.45) 



where ; is the vector of second-order operators de- 
fined in Eq. (|4.42[) . It is easily shown that the first term 
cancels with the contributions coming from the second- 
order operators, and that the contributions from second- 
order wave function are exactly zero. Hence, that block 
is simply 



,(2) 
1,1 



(0(°)|[Tl,[iT«,Tt]|0«) 



(4.46) 



which makes it frequency- independent. Its calculation 



gives 



r r (2)i 

L 1 1 i\kc,ia 



. M kd M di ^ M la M lc 



(4.47) 



Sacsr^ (le\\kd)(dl\\ei) 

Ide 



9. 2^ 



^im,de 

(ld\\mc)(dl\\ma) 

^Im.ad 



Imd 

M ak M ld , M ci M ka 



M d k{ad\\ci) 



M lc (lk\\ai) 



I 

y-^ (ce\\ad)(di\\em) 

, ^im.de 
ma 

y-r (ce||mz)(aA;||me) 

me 



-y 



ml 



(ce\\ad)(dk\\ei) 



£ik : de 



(ik\\ml){ac\\ml) 



(4.48) 



The block T\,2 and its adjoint is of at least first-order, 
due to the fact that the space is orthonormal. For that 
reason, it is not affected by the orthonormalization at 
this level of approximation. Its calculation gives 

2) 



2,Ufec,jfei a = ~ s ik(bc\Wj) + S jk (bc\\ai) 
- 5 bc {ai\\kj) + 5 ac (bi\\kj) 







(4.49) 



2,lJ ck t jbia 

Finally, the block T2.2(w) gives 

( w )]ldkc,jbia = (uj — e-ij >a b)5jl5ik&ca.bdb 

[r 2 2 !(o;)] cfcrf ;, i6ia = (4.50) 

Notice that double excitations are treated only to zeroth- 
order in a second-order approach. To obtain a consistent 
theory with first-order corrections to double excitations, 
one should go at least to third order. This however be- 
comes computationally quite heavy. 
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It is interesting to speculate what would happen if we 
were to include the first-order doubles correction within 
the present second-order theory. There are, in fact, indi- 
cations that this can lead to improved agreement between 
calculated and experimental double excitations, though 
the quality of the single excitations is simultaneously de- 
creased due to an unbalanced treatment. (63l. I&H 

We can now construct the PP necessary to construct 
the 2nd approximation of the xc-kernel Eq. (|2.46p accord- 
ing to Eq. (|4.26p . Since the the localizers of both left- 
and right-side are constructed from the non-interacting 
KS PP, we are only concerned with ph and hp contri- 
butions. This means that the blocks involving pp or hh 
indices, corresponding to density shift operators, can be 
ignored at this level of approximation. This simplifies the 
construction of P(w) in Eq. (|4.26[) , which up to second- 
order gives 

n (o+i+2),-i (w) = (T 1 t|T 1 t )- 1 p(°+ 1 + 2 )( W )(T 1 t |T 1 + )- 1 . 

(4.51) 

Separating ph and hp contributions, the PP takes the 
form of a 2 x 2 block-matrix in the same spirit as the 
LR-TDDFT formulation of Casida, 



n (0+l+2),-l (w) = 

(\ o \ / p(°+ 1+2 )H p(°+ 1+2 )H \ / i o 
l,o -l J ^ p( + 1 + 2 )( w ) p(°+ 1+2 )H )\o -l 

/ p(0+l+2) (w) _p(0+l+2)( w ) \ 
- ^ _p(0+l+2)( w ) P(°+1+2)( W ) J ■ 

It follows that, 

n; 1 H-n( 0+1 + 2 )- 1 H = 

/pd+^H -if+ 2) \ 
I _ r a+2) p ( i+2) (w) I • 



(4.52) 



(4.53) 



Note that the off-diagonal (ph,hp)- and (hp,ph)-blocks 
are frequency-independent and that the diagonal blocks 
are given by Eq. (|4.44j) . Ignoring localization for the 
moment, we may now cast the present Kohn-Sham 
based second-order polarization propagator approxima- 
tion (SOPPA/KS) into the familiar form of Eq. ([Pj) 
with, 



A iatjb (u) = 8 itj S atb (e a - ei) + P&ffiw) 



Bi atbJ (uj) - - (r^ 2 - 1 



ia.bj 



(4.54) 



by mixing the P^ 1+2 \uj) and T^ 2 ^ terms, 

Ai a jb{w) = 5i,jSa,b\£a — e i) 

+ [(A s ) ftp ,^( W )p( 1+2 )( W )(At), iPi ^( w ) 
+ [(A s ), ip ^( w )p( 1+2 )( W )(At) p ,^ p ( w ) 
(A s ) hp , ph (w)r( 1+2 \At) hp , hp (u) 
(A s ), ip , hp ( w )r( 1+2 )(At) ph ^ p ( w ) 
B ia ,bj(u) = (A 3 ) hp>hp p( 1+2) {uj)(Al) hp>ph 
+ [(A a ) hPiPh P< 1+a >(w)(At) phiI , h 
(A s ) hp , ph (w)r( 1+2 \At) hp , ph (u) 

■ "(A s ), iPihp ( w )r( 1+2 )(At) p ^ p/i ( w ) 



ia,jb 
ia,jb 



ia,jb 
ia,jb 
ia,bj 
ia,bj 
ia,bj 
ia,bj 



(4.55) 



Of course this extra complication is unnecessary if all 
we want to do is to calculate improved excitation en- 
ergies and transition amplitudes by doing DFT-based 
many-body perturbation theory. It is only needed when 
our goal is to study the effect of localization on purely 
TDDFT quantities such as the xc-kernel and the TDDFT 
vectors X and Y. 



D. Tamm-Dancoff Approximation 

The present theory, without the localizer, takes a 
particularly transparent and computationally convenient 
form when the Tamm-Dancoff approximation (TDA) is 
made. This SOPPA/KS + TDA consists of consider- 
ing only the (hp,hp)-block. The excitation energies and 
transition amplitudes are found by solving the pseudo- 
eigenvalue equation, 



1,2 



2.2 



L (l) 

2,1 



or the equivalent true eigenvalue equation, 



Localization [Eq. (|2.46|) ] will complicate these formulae 



4(0+1+2) 4 (i) 
""1,1 -^1,2 



A {1) 
^2,1 



,(0) 



C 2 



Co 



Ci = wCi , 
(4.56) 

(4.57) 
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Here 



This separates into the static part, 



(< +1+2) ) t = S, k F^ - *a,cif fe +1+2) + (0i||*c) 
\ / kc,ia 

( A 2i) = -Si,k{bc\\aj) + 5 jtk (bc\\ai) 

\ ' / kc,jbia 

- S b , c {ai\\kj) + 5 k ,j(U\\kj) 

Si,kS c ,aSd,be a b,ij , (4.58) 



(43) 



ldkc,jbia 



where Fr°s +1 ^ = 5 rtS e r + M^ c s is the matrix of the Hartree- 
Fock operator constructed with Kohn-Sham orbitals and 

M La M, c 



1 ^-v (ld\\mc)(dl\\am) 
M ktd M d A 



,(o+i+2) = ^,(0+1) 

^k % .k j^^^j 



1 

2 . J ^im,de 



£i,d 

1 ^ (Je||fcd)(dJ||ei) 

Ld,e 



(4.59) 



include second-order corrections. (Note that extra fac- 
tors of 1/2 will occur in these expressions when spin is 
taken explicitly into account.) It is immediately seen 
that truncating to first order recovers the usual CI sin- 
gles (CIS) equations in a noncanonical basis set. 



V. DRESSED TDDFT 

Dressed TDDFT (D-TDDFT) is a way to mix conven- 
tional adiabatic TDDFT with a nonadiabatic correction 
derived from many-body theory. Now that we have all 
the pieces, let us assemble a first-principles PP version of 
D-TDDFT. 

The first step along the way is to understand how 
to avoid double counting of correlation effects. That is 
we must generalize the Maitra, Zhang, Cave, and Burke 
(MZCB) ansatz to more than one doubly-excited state. 
A perturbative generalization of the formula (|1.23[) in- 
volving a single doubly-excited state to the situation in- 
volving a manifold of doubly-excited states is, 



UJ S 



E 

D' 



D' 



0J — LOjji 



(5.1) 



\ X 



D 



LV 

1 Ujy 



and a remaining dynamic part, 



una 



uj — ujd' 



(5.2) 



(5.3) 



Since adiabatic TDDFT must agree with static DFT re- 
sponse theory, we can easily identify ojaa &s the part 
described by conventional adiabatic TDDFT. 

In order to treat the remaining part, we restrict our- 
selves to the situation that the excitation energies are 
well separated with a particular w ujg . In this case, 
we can write that, 



una 



d 



ujd I u>- uj0 



■n 



oj — ojd 



(5.4) 



This later formula is the dressing to be obtained from 
PP theory. Note that it assumes that the true excitation 
energy, u » ujd. 

Equation (|4.56[) is of just the right form to apply this 
type of ansatz. Assuming just one (or at most a few) dou- 
ble excitations, our MZCB-like PP ansatz is the pseudo- 
eigenvalue equation, 



with 



1-1,2 



(wl 



2,2 



2.1 



C x = uC x , (5.5) 



kc,jbia 



0$) 



ldkc,jbia 



6i,kSa,c£a,i + (M/ffxcl*°) 

-S l . k (bc\\aj) + 5j t k(bc\\ai) 
$bA ai \\ k j) + h,j(fii\\kj) 

Si,kS c ,aSd,b^ab,ij ■ 



(5.6) 



Equation (|5.5[) is equivalent to the computationally more 
convenient form, 



A M Ai )2 

A-2,1 A2.2 



Co 



Ci 

Co 



(5.7) 



which is what we use in practice. 

The restriction to one (or at most a few energet- 
ically similar) double excitations can be justified on 
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two grounds. First, adiabatic TDDFT already contains 
electron correlation effects and we wish to avoid dou- 
ble counting. Second, the Gonze-Schefffer relation tells 
us that the localizer transforms diagonal ia-matrix el- 
ements of the spatially-local kernel, fH XC (l,2;u), into 
matrix elements of the spatially nonlocal operator, 
Kh xc (1, l';2,2';w), at frequencies, u = e a ,i- Were we to 
consider an increasingly large number of 2h2p-excitations 
in our theory, there may come a point where it will be 
necessary to include the localizer when interfacing PP 
and TDDFT formulae. For the moment we do not do 
this in any of the calculations presented in this paper. 
It may be worth emphasizing that the localizer is un- 
necessary when carrying out PP based upon the Kohn- 
Sham orbital hamiltonian as zero-order estimate, but is 
expected to become necessary at some point when inter- 
facing the spatially-local TDDFT and spatially nonlocal 
PP formalisms. 

The formulae (|5.6p are highly reminiscent of Casida's 
earlier work on a PP correction to adiabatic TDDFT. 
These formulae are just the TDA form of the order- 
consistent development of Eq. (3.9) in Ref. (22|. Com- 
paring with Eqs. (|4.58[) and (|4.59[) , we see that the 
block of lhlp-excitations is given by conventional adi- 
abatic TDDFT while the missing terms involving 2h2p- 
excitations are taken over from SOPPA/KS+TDA the- 
ory. We refer to them here by the acronym D-TDDFT 
TDA. It is troubling that the zero-order estimate of the 
double excitation energy is just given by KS orbital en- 
ergy differences in this order-consistent theory. These 
may or may not be accurate enough in practice, so we also 
introduce the notion of extended dressed (ED) TDDFT 
in which, 

(Ag£) = {9 K s\ckWMjaU\9 KS ) 

V / ldkc,jbia 

- {^ks\H\^ks)- (5.8) 

An explicit expression is straightforward to work out with 
Wick's theorem but is long and so will not be given here 
[see Eqs. (49) and (50) of Ref. [TO]]- 

Let us now turn to a specific example, namely that of 
butadiene (CH2=CH-CH=CH2), whose excitation spec- 
tra have been extensively studied in the past (see e.g., 
Ref. (6f| and references therein). Let us examine the 
electronic structure of this molecule in a qualitative man- 
ner before embarking on a quantitative study. Fig- 
ure IVl shows the frontier molecular orbitals (MOs) of this 
molecule, easily obtained by qualitative arguments famil- 
iar to most chemists. Several possible excitations may 
be constructed from these MOs with the dominant ones 



FIG. 5. Frontier molecular orbitals of butadiene calculated 
using the local density approximation. 




shown in Fig. |V1 We can anticipate a fairly pure one- 
hole one-particle (lhlp)-excitation of 1 1 B U symmetry. 
However excitations of 1 A g symmetry are more com- 
plicated, since it is possible to construct two energetic 
nearly degenerate lhlp-excitations from the HOMO — >• 
LUMO+1 and HOMO-1 -> LUMO transitions. (HOMO 
and LUMO are standard acroynms respectively for "high- 
est occupied molecular orbital" and "lowest unoccupied 
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FIG. 6. Dominant contributions to the low-lying excited 
states in trans-butadiene. Huix-Rotllant and Casida: 
dressed TDDFT 



TABLE I. Butadiene vertical excitation energies (eV). 



lBu 



2 b. 



2A< 



2a u 



lb. 



A A 



lau -LI 



molecular orbital.") These will interact to form symmet- 
ric and antisymmetric combinations, one of which may 
strongly mix with the HOMO -> LUMO double (2h2p) 
excitation. The bright 1 l B u state is experimentally rel- 
atively easy to locate with absorption spectroscopy while 
the location of the spectroscopically dark 2 1 A g is much 
harder to locate experimentally. Recent experiments 
seem to indicate that the 2 1 A a state is slightly lower in 
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while theory is divided 
state is slightly higher 66] or 



energy than the 1 1 B U state, 
as to whether the 2 1 A g 

slightly lower [65| than the 1 1 B U state (Table HJ). Prob- 
ably at this point we should only conclude that the two 
states are energetically quasidegenerate. 

Butadiene has been treated at least twice in the past 
using D-TDDFT. [H, [H The earliest treatment was by 
Cave, Zhang, Maitra, and Burke [l9j]. They first carried 
out TDDFT calculations to find the singles-excitation 
energy us, and the corresponding single-excitation con- 
figurations. They then added an extended dressed term 





2M 9 


1 B u 


Present Work a 




TDLDA TDA 


6.42 


coo 


D-TDLDA TDA 


3.91 


6.38 


ED-TDLDA TDA 


6.02 


6.38 


10-ED-TDLDA TDA 


6.253 


6.249 


Cave el al 


b 




TDB3LYP 


".02 




TDPBEO 


7.22 


5.98 


D-TDPBEO 


5.93 




D-TDPBEO TDA 


6.28 




Mazur and Wlodarczyk 


c 


TDB3LYP TDA 


6.56 


t\ no 


TDPBEO TDA 


6.77 


6.13 


D-TDB3LYP TDA 


6.06 


6.02 


D-TDPBEO TDA 


6.34 


6.13 


Reference Values 




SOPPA^ 


7.49 


6.00 


ADC(2)-s~^ 


8.34 


7.09 


ADC(2)-x3 


5.19 


6.69 


CASPT2 - £" 


6.27 


6.23 


Best Estimate^ 


6.55 


6.18 


Experiment^ 


5.4-5.8 


5.91 



a Present work: 6-311G(d,p) basis. 

b Calculations of Cave, Zhang, Maitra, and Burke from Ref. [Ts| : 

6-311G(d,p) basis. 
c Calculations of Mazur and Wlodarczyk from Ref. [Soli : 

6-311++G** basis. 
d Strict (s) and extended (x) ADC(2)/6-31G calculations from 

Table 2 of Ref. [H. 
° From Ref. 
f From Ref. 

s Values taken from Ref. |65|| . See also the discussion of these 
experimental values in that reference. 



into both the A and B matrices and solved the full 
and TDA dressed TDDFT problem in a noniterative 
perturbation-like manner using electron repulsion inte- 
grals constructed with Hartree-Fock orbitals. Note that 
their double-excitation energy , ojd, is calculated not in 
terms of orbital energy differences but in terms of the 
difference between the Hartree-Fock ground state and 
the energy of a doubly-excited Hartree-Fock determi- 
nant, similar to our ED-TDDFT theory. Their results are 
shown in Table Q] While the TDPBEO 1 1 B U excitation 
energy is reasonable, the TDPBEO 2 1 A g excitation en- 
ergy is much too high compared with the 1 1 B U excitation 
energy. Their dressed TDDFT gives &2 1 A g energy simi- 
lar to the 1 1 B U excitation energy. Mazur and Wloarczyk 
[2pj ] have recently improved on the work of Cave et al. by 
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calculating their electron repulsion integrals from Kohn- 
Sham, rather than Hartree-Fock, orbitals and by car- 
rying out an iterative solution of an improved form of 
Eq. (|1.23|) . They obtained results remarkably similar to 
those of Cave et al. Note that none of these dressed calcu- 
lations affects the energy of the \ X B U state which belongs 
to a different irreducible representation of the molecular 
point group. 

We have carried out D-TDDFT calculations with our 
PP formalism using the Grenoble development version 
of the DeMon2k program [69|. The implementation of 
TDDFT in deMon2k has already been described (70} . 
The implementation of D-TDDFT and extensive tests 
of this option will be described elsewhere. [7l| Suffice it 
to say that calculations were carried out at the experi- 
mental geometry using the same orbital basis set as in 
used by Cave, Zhang, Maitra, and Burke 0]. Our cal- 
culations are limited to the local density approximation 
(LDA) rather than the hybrid functionals (PBEO and 
B3LYP) used in the previously mentioned work. Our re- 
sults are collected in Table U and compared against other 
D-TDDFT calculations and benchmark values from the 
literature. We confirm the previous result of Hsu, Hi- 
rata, and Head-Gordon [72[ that non-hybrid function- 
als give reasonable agreement with benchmark values for 
the 2 1 A g excitation energy. As seen in the table, hy- 
brid functionals are even worse in this respect. Including 
the lb 2 n 2al double excitation at the D-TDLDA TDA 
level of theory lowers the energy of the 2 1 A g state by 
a very large 2.5 eV. This is because the simple orbital 
energy difference of the 2h2p excitation energy at 8.01 
cV is much too close to the 2 1 A g excitation energy (6.42 
eV) obtained without the doubles correction. Using the 
extended doubles formula (|5.8p with a 2h2p excitation 
energy of 23.84 eV leads to a much more reasonable low- 
ering of the TDLDA TDA energy by 0.4 eV. As explained 
above, the energy of the 1 1 B U state remains unchanged 
by the introduction of a double excitation of l A g sym- 
metry. 

It should be emphasized that our implementation of 
ED-TDDFT is a still more general implementation of the 
previous implementation of dressed TDDFT by Mazur 
and Wloarczyk [20} which was already an improvement 
on the work of Cave et al. [l9[ Not only do we obtain 
the equivalent of the iterative answer with a noniterative 
algorithm [compare Eqs. (|5.5[) and ()5.7[) ] but we are able 
to include a few additional double excitations. (More- 
over we have given the formal framework for further sys- 
tematic improvements of our theory.) Thus we thought 
it interesting to go a little further in this preliminary 



study by including the hrst 10 doubly-excited states in 
the ED theory (without including the localizer). This 
is the calculation denoted 10-ED-TDLDA TDA in Ta- 
ble HI This introduces some double excitations of 1 B U 
symmetry. While some double counting of correlation is 
inevitable, 10 double excitations is still a very small num- 
ber compared to the number usually used to describe dy- 
namic correlation. Both the 2 1 A g and 1 1 B U states have 
shifted to nearly the same value, in better agreement with 
reference values. 

Since our version of D-TDDFT is based upon SOPPA 
and strict ADC(2) formulae, it may be worth compar- 
ing against these values. Table Q] shows the results of a 
SOPPA calculation that we obtained with the Dalton 
code [73} using the same basis and geometries as our D- 
TDDFT calculations. The 1 1 B U SOPPA excitation en- 
ergy is similar to that obtained with the TDLDA TDA 
and D-TDLDA TDA. However, unlike in our D-TDLDA 
TDA calculation, the SOPPA 2 1 A g excitation energy re- 
mains higher in energy than the 1 l B u excitation energy. 
This is consistent with the much larger double excita- 
tion orbital energy difference when Hartree-Fock, rather 
than Kohn-Sham, orbital energies are used. We also 
give strict JADC(2)-s] and extended [ADC(2)-x] ADC(2) 
from Ref. j67|. The ADC(2)-s result should be compara- 
ble to our SOPPA result but overestimates our result 
for the excitation energies of both states because the 
ADC(2)-s result was calculated with a less extensive basis 
set. The ADC(2)-x method is to the ADC(2)-s method 
what our ED-TDLDA TDA method is to our D-TDLDA 
TDA method. Interestingly, the ED-TDLDA TDA 2 1 A g 
energy is greater than the D-TDLDA TDA 2 1 A g en- 
ergy while the ADC(2)-x 2 1 A g energy is smaller than 
the ADC(2)-s excitation energy. This is consistent with 
the correction to the simple orbital energy diffrencc esti- 
mate being very different for the ED-TDLDA TDA and 
ADC(2)-x schemes. In both cases the extended method 
brings the calculated results closer to experiment. 

In finishing this section, it is best to recall that the re- 
sults given here are only indicative. A much more thor- 
ough study will be reported elsewhere. 71] In particular, 
it should be kept in mind that the 2 1 A g and 1 1 B U states 
lie above the TDLDA ionization threshold which is artifi- 
cially low (5.92 eV). [74| Nevertheless the present results 
give a first indication that our method is up and running 
and ready for more extensive testing. A particularly in- 
teresting question will be to what extent our mixture of 
TDDFT and PP methods will break symmetry-imposed 
state degeneracies. 
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VI. CONCLUSION 

This article has presented a detailed discussion of 
how the frequency-dependent spatially-local xc-kernel of 
TDDFT may, in principle, be obtained from the corre- 
sponding spatially-nonlocal quantity in PP theory in such 
a way as to take full account of the interaction between 
any number of single- and higher-order excitations. This 
may be termed ab initio TDDFT. 24| Because we have 
made extensive connections with other similar work in 
the literature, it seems best to conclude by emphasizing 
what is new in the present work. 

The first thing that is new is that we have clari- 
fied the relation between the PP work of Casida [22[ 
and similar work within the solid state physics commu- 
nity based upon the Bethe-Salpeter equation (BSE) (33 - 
l39l |42| . This is important because different approxi- 
mations have been developed within the PP approach 
in chemical physics for molecules than those developed 
within the BSE approach for solids, yet there is a great 
deal of interest in these two communities in sharing infor- 
mation with the aim of eventually developing a physical 
approximations applicable, say, at the nanoscale interface 
between molecules and solids. 

The second thing that is new in the present work is that 
we have corrected the previous PP work of Casida [22 1 
to include the effect of spatial localization. Within the 
Born approximation, the resultant expression for the xc- 
kernel factors into a product of dynamic localizers and 
the nonlocal PP analogue of the xc-kernel. This localizer 
arises naturally in OEP and TD-OEP methods. Inter- 
estingly the dynamic localizer represents the fact that 
spatial localization necessarily implies an additional fre- 
quency dependence (delocalization in time). Near a pole, 
when the excitations are well-separated (which is nor- 
mally what happens for molecules), we obtain a par- 
ticuarly easy demonstration of the Gonzc-Schcfflcr rela- 
tion 5l| which is, among other things, a key ingredient 



in a correct treatment of charge-transfer excitations in 
TDDFT. [52j With a future implementation in mind and 
explicit exploration of the properties of the localizer, we 
also briefly discussed finite-basis representations of the 
localizer in an appendix. 

The third thing that is new here is that we have de- 
veloped explicit formulae similar to SOPPA or ADC (2) 
formulae but based upon a Kohn-Sham zero-order hamil- 
tonian. This involves resumming more than one hundred 
correction terms but results in compact formulae in the 
end. We have shown that previous EXX formulae re- 
ported in the literature are confirmed exactly by trun- 



cating our formulae to lower order. 

The fourth thing that is new is that we have discussed 
and given what we believe to be an improved formulation 
of dressed TDDFT. A preliminary calculation on buta- 
diene confirms that it is indeed necessary to go beyond 
a zero-order description of double excitations to at least 
a first-order description. The results confirm that the 
method is correctly implemented and are sufficiently en- 
couraging that we are in the process of carrying out a 
benchmark study of this approach on a much more ex- 
tensive set of molecules. [7l| 

Finally we have presented all the pieces in this paper to 
go from ab initio PP theory to extract the corresponding 
frequency-dependent xc-kcrncl of TDDFT. The formal- 
ism itself offers insight into the analytic structure of this 
xc-kernel and why dressed TDDFT works and eventually 
why it might fail and could be further improved. Actual 
calculations of the ab initio xc-kernel can be expected 
to lead to additional insight, but the necessary code will 
take some time to put into place. 
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Appendix A: Collapse, Expansion, and Localization 
Operators 

A problem which shares many aspects in common with 
the present work is that of finding the optimized effective 
potential (OEP). HKzl] The OEP is the local potential 
whose orbitals minimize the Hartree-Fock energy expres- 
sion, and is, up to a linear response approximation, iden- 
tical to the exact exchange-only Kohn-Sham potential, 
V x , whose functional derivative is the exact exchange- 
only kernel, f x . In using Harriman's collapse operator 
in the present work, we have chosen to follow Hamel, 
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Casida, and Salahub who pointed out that this opera- 
tor allows a p articularly elegant formulation of the OEP 
problem. |53j The notion of a collapse operator emerges 
naturally in this context. This appendix introduces the 
notion of collapse, expansion, and collapse operators in 
the static case. Some exact and approximate properties 
are also reviewed. 

Harriman considered the relation between th e sp ace of 
kernels of operators and the space of functions. [Zs[73] In 
order to maintain consistency with the rest of this paper, 
we will generalize Harriman's notion from space-only to 
space and spin coordinates. Then the collapse operator 
is defined by, 

TS X (1 I 2) = E 1B (1,1). (Al) 

The adjoint of the collapse operator is the so-called ex- 
pansion operator, 

ft« x (l)=« ai (l)<y(l-2). (A2) 

The notions are general even if the above two equations 
suggest that we apply the collapse operator to the kernel 
of the exchange-only self-energy (i.e., the Hartree-Fock 
exchange operator) and that the expansion operator be 
applied to the exchange-only potential, which is our in- 
tention. 

The exchange-only OEP problem may be written as, 

(h H + v x ) ^ EP (x) = e? EP ^ EP (x) , (A3) 

which resembles the Hartree-Fock problem, 

(h H + E a ) (x) = ef F Vf F (x) ■ (A4) 

It turns out that the answer to the OEP problem is to 
require that the potential v x be chosen so that the pertur- 
bation transforming the OEP problem into the Hartree- 
Fock pro blem makes the linear response of the density 
zero. [25l . |75| Since the definition of the exchange-only 
Kohn-Sham potential is that it produce the Hartree-Fock 
density, it is easy to see that the OEP v x is simply a linear 
response approximation to the exact Kohn-Sham v x . 

We can express this linear response condition very el- 
egantly using the Harriman operators as, 

= in s (E a - T+t;,) . (A5) 

In order to keep the notation compact by suppressing in- 
tegrations, we have introduced a matrix notation: The 
operators T and the OEP (Kohn-Sham) independent 



particle response function (retarded form of the polariza- 
tion propagator), II S , are represented as matrices. The 
operators H x and v x are represented as column vectors. 
Thus Eq. (|A5|) becomes, 

ru s r^ Vx = Tn s s x , (A6) 

which may be recognized as a weighted least square fit of 
T^v x to Ti x . It may be solved, 

v x = A S T, X , (A7) 

where, 

A s = (Tn s Tt) _1 Tn s , (A8) 

is the localizer. (In particular it is the single-particle 
localizer because it depends upon the single particle re- 
sponse function II S .) Note that this only makes sense if 
care is taken to eliminate the null space of YII S before 
taking the inverse (i.e., singular value decomposition.) 

One of the strong points of this formulation of the OEP 
problem is that it lends itself to a spectral space repre- 
sentation in terms of orbital and auxiliary basis sets. [53[ 
However solving the OEP equation in this manner is not 
without its pitfalls. Hamel, Casida, and Salahub empha- 
sized that, 

"The v x which emerges from finite basis OEP 
calculations is not strictly unique, though its 
matrix is uniquely defined up to an additive 
multiple of the matrix of the identity." [53j j 

In fact, an improper balance between orbital and aux- 
iliary basis sets can lead to unphysical wiggles in v x . 
Staroverov, Scuseria, and Davidson have pointed out that 
an improper balance between the orbital and auxiliary 
basis sets can also lead to Hartree-Fock e nerg ies which 
are incorrectly below the true OEP energy. [78[ Thus the 
orbital and auxiliary basis sets should be regarded not as 
independent but rather are coupled. Methods for avoid- 
ing this pitfall have been discussed in the recent litera- 
ture [79l - [81| and some are described in the following para- 
graphs. However, for the purposes of the present work, 
we need only assume that the pitfall may be avoided us- 
ing appropriately balanced basis sets. 

One way to understand the localizer is to examine its 
static limit. There are several related approximate so- 
lutions to the static OEP equation [Eq. (|A6[) ]. We fo- 
cus here on the Krieger-Li-Iafrate (KLI) approximation 
because of its physical transparency. However a draw- 
back of the KLI approximation is lack of unitary in- 
variancc which may be removed by a suitable general- 
ization giving what is variously known as the common 
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energy denominator approximation (CEDA) [82[, local- 
ized Hartree-Fock (LHF) [83| , or effective local potential 
(ELP) method [zl, H(|, depending upon who derived it. 

Krieger, Li, and Iafrate [HI, HH followed up on a foot- 
note in Sharp and Horton's seminal OEP paper. [75[ 
Amazingly the resultant KLI expression works almost as 
well as the original OEP method. The trick is to recog- 
nize that the generalized susceptibility may be rewritten 
as, 



II S (1, 2; 3, 4) = ]T]T 

i p^i 
occ 

+ EE 



&(lW;(2)#(3ty„(4) 



Vv(l)^*(2)</>;(3)^(4) 



l -,P 



(A9) 



because the terms where p corresponds to an occupied 
orbital all cancel each other. Now make use Unsold's 
approximation j86j in which we replace the excitation 
energies with a single average excitation energy, 



n s (l,2;3,4) = -^^ 

i p^i 



^(l)V*(2)V?(3)Vp(4) 



EE 

i Py^i 



Ae 

^(1)^(2)^(3)^(4) 
Ae 



(A10) 

(Csizmadia and Sylvain gave evidence to suggest that 



J_ 1 ™ 1 

Ae _ N nr .r, ^ U. ' 



(All) 



is a good choice when applying Unsold's approximation 
to the calculation of polarizabilities H3| ■) Using the com- 
pleteness relation, 



£>p(1)W)=*(1-2), 



(A12) 



then 



11.(1, 2; 3, 4) = -— ( 7 (1, 3)5(4 - 2) + 7 (4, 2)<5(1 - 3) 
Ae 



2^^(1)^(2)^(3)^(4) 



(A13) 



Inserting into Eq. (|A6|) and assuming real orbitals yields 
after a little algebra, 



KLI 



(1) 



piX) 



p(i) 



(A14) 



Notice how Ae has completely canceled out! The first 
of the two terms is the Slater potential, [1^ while the 
second of the two terms is a "derivative discontinuity 
correction." 

The presence of the exchange potential on the right 
hand side of Eq. (|A14[) shows the real complexity under- 
lying the apparently simple equation (|A7[) . Nevertheless 
we have here an implicit approximate expression for the 
static localizer. A similar statement also holds for the 
CEDA (LHF, ELP) approximation. 



Appendix B: Finite-Basis Representation of the 
Dynamic Localization Operator 

We discuss the problem of carrying out dynamic local- 
ization, 



v(l,u) = A(w)£(l,2;w) 



(Bl) 



[e.g., Eq. p.6p ] in a finite basis representation. Such a 
computational approach is probably the most rigorous 
direct way to explore properties of the dynamic localizer. 

In a typical finite basis method, v is expanded in a 
finite auxiliary basis set, 



«(l;w)=5^ fl i(l)6/(a;). 



(B2) 



It is evident that the gi must be completely contained 
within the space of orbital products, otherwise erratic re- 
sults will be obtained when expanding the orbital space. 
The implication is that, at a minimum, there is a dou- 
ble convergence problem. For each auxiliary basis set, it 
is necessary to expand the orbital basis set until conver- 
gence is reached. Then the auxiliary basis set can be ex- 
panded and convergence attained once again with respect 
to the orbital basis set. The reverse order (convergence 
of the auxiliary basis set before the orbital basis set) is 
not possible. This is at the heart of problems mentioned 
in Appendix HI (zUll 
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Nevertheless, with these difficulties in mind, let us 
make the expansion in the auxiliary basis set [Eq. (|B2|) ] 
and rewrite the localization equation as, 



J2x I , J {oj)bj(uj) = Y I (u 



(B3) 



with 



i a 
occ virt 

-EE 



(lpa\gi\lpi)(lpi\gj\lpa) 



{^i\gi\^a){^a\9j\ll>i) 



(B4) 



and 



w = EE 

i a 
occ virt 

-EE 



(lpa\gi\lpt)('lpi\^(L0)\lpa) 
(lpi\gi\lpa){'lpa\^(uj)\lpi) 



(B5) 



[Here we are using the compact notation introduced in 
Eq. (jTTTJ) .] For real basis sets, 



2r„ 



?. a 
occ virt / 



uj 2 - e 2 ■ 



(^a\gi\ipt)(A\gj\ipa) 



w=EE 



2e„ 



^ 2 -e 2 



{lka\9l\i>i){^i\^{^)\^a) 



(B6) 



The solution to the localization problem is, 

,7 

1,7 

= 5^ S j(i)A Jiia (w)(Vi|x:(w)|Va) 



(B7) 



with the localizer supermatrix given by, 

2e a ,i{ipa\gi\ipi) 



a. 1 m = y^xj}(u) r- 



uj 2 -e 2 ■ 

a.i 



(B8) 



Similar equations can be written down for the KLI 
and CEDA/LHF approximations but the critical prod- 
uct space is then that of products of occupied orbitals, 
i>iipj, which seems to be somewhat easier to handle, to 
judge by the success of the ELP method [zl,[8^. For ex- 
ample, if we make the CEDA/LHF approximation with 
the energy dependent n s (w), then there will be a factor 
of 2Ae/[w 2 — (Ae) 2 ] which will cancel out and we will 
arrive back at the CEDA/LHF equations but this time 
keeping the frequency dependence in v x (u) and in 



w(l;w) 



p(i) 

S°,j C V'i(i)(V'iK^)-sMIV'j)^(i) 



(B9) 

(This is a more general result than the KLI approxima- 
tion described in Appendix |A] The KLI approximation 
is obtained when i = j in the second sum.) The matrix 
equations are, 

Y,X I ,jbj{w)=Y J {Lo) 
j 

occ occ 

Xi,J = ^2(ipi\9i9j\ipi) - ^{^i\9i\^j){^j\9j\^i) 

i i-.j 

occ occ 

y t (u,) = J2(i>i\9i£(uM) -^<V>i|<?i|^><V#MIV'i> 

i i,j 

bj{w) = Y,Xj}Y I {w). (BIO) 
I 

It is to be emphasized that the localizer is frequency in- 
dependent in this approximation. Although this may 
seem like a great computational simplification, it is at 
the price that the Gonze-Scheffier relations [Eq. (|3.14[) 
and Eq. ()3.15[) ] are no longer exactly satisfied. 

A frequency-dependent CEDA has recentl y be en intro- 
duced recently by Gritsenko and Baerends. [2l[ 

Appendix C: Diagrammatic Rules 

The representation of MBPT expansions in terms of 
diagrams is very convenient for bookkeeping purposes. 
Indeed certain ideas such as the linked-cluster theorem 
fo2j or the concept of a ladder approximation (see e.g., 
Ref. [46| p. 136) are most naturally expressed in terms 
of diagrams. Perhaps most importantly diagrams drawn 
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,{t,t') = 9(t 




where, 



FIG. 7. Basic time-ordered finite basis set representation PP 
diagram. 



according to systematic rules allow an easy way to check 
algebraic expressions. Several types of MBPT diagrams 
exist in the literature. These divide into four main 
classes which we call Feynman, Abrikosov, Goldstone, 
and Hugcnholtz. Such diagrams can be distinguished by 
whether they arc time-ordered (Goldstone and Hugcn- 
holtz) or not (Feynman and Abrikosov) and by whether 
they treat the electron repulsion interaction as a wavy 
or dotted line with an incoming and an outgoing ar- 
row at each end (Feynman and Goldstone) or in a sym- 
metrized way as a point with two incoming and two out- 
going arrows (Abrikosov and Hugenholtz). These dif- 
ferences affect how they are to be translated into alge- 
braic expressions as does the nature of the quantity be- 
ing expanded (wave function, one-electron Green's func- 
tion, self-energy, polarization propagator, etc.) Given 
this plethora of types of diagrams and the difficulty of 
finding a clear explanation of how to read polarization 
propagator diagrams, we have chosen to present an ap- 
pendix summarizing how our diagrams should be trans- 
lated into algebraic expressions. This is perhaps espe- 
cially necessary because while the usual practice in the 
solid-state literature is to use time-unordered diagrams 
with electron repulsions represented as wavy or dotted 
lines (i.e., Feynman diagrams), while the usual practice in 
the quantum chemistry literature of using time-ordered 
diagrams with electron repulsions represented as points 
(i.e., Hugenholtz diagrams). Since we first derive alge- 
braic expressions and then write the corresponding dia- 
grams, rather than vice versa, we will content ourselves 
in this appendix on how to translate diagrams into al- 
gebraic expressions, rather than upon how to generate 
complete non-redundant sets of diagrams. We will de- 
velop the rules gradually, including examples along the 
way. 

The PP expressed in an orbital basis is, 
11(1,2,3,4;*-*') = ^ n sr ,, p (t-t')Vr (1)^(2)^(3)^(4). 

pqrs 

(Cl) 



-iO(t - t')(0\f^(t)s H (t)q^(t')p H (t')\0) 

l9(t' - t)(0\ql(t')p H (t')f^(t)s H (t)\0) . 

(C2) 



This makes it clear that the PP is a two time particle- 
hole propagator which either propagates forward in time 
or backward in time. To represent it we introduce the 
following rules: 

(1) Time increases vertically from bottom to top. This 
is in contrast to a common convention in the solid- 
state literature where time increases horizontally 
from right to left. 

(2) A PP is a two time quantity. Each of these 
two times is indicated by a horizontal dotted line. 
This is one type of "event" (representing the cre- 
ation/destruction of an excitation). 

(3) Time-ordered diagrams use directed lines (arrows) . 
Down-going arrows correspond to holes running 
backward in time, that is, to occupied orbitals. Up- 
going arrows correspond to particles running for- 
ward in time, that is, unoccupied orbitals. 

At this point, the PP diagrams look something like Fig. [7] 
Fourier transforming leads us to the representation shown 
in Fig. IS] An additional rule has been introduced: 

(4) A downward u> arrow on the left indicates forward 
ph-propagation. An upward uj arrow on the right 
indicates downward ph-propagation. 

Diagrams for the corresponding position space represen- 
tation are shown in Fig. |H1 Usually the labels (p, q, r, 
and s or 1, 2, 3, and 4) are suppressed. If the ui arrows 
are also suppressed then there is no information about 
time-ordering and both diagrams may be then written as 
a single time-unordered diagram as in Fig. 1101 Typical 
Feynman diagrams are unordered in time. 

Perturbation theory introduces certain denominators 
in the algebraic expressions corresponding to the dia- 
grams. These may be represented as cuts between events. 

(5) Each horizontal cut between events contributes a 
a factor (±ui + *£ p e P - Eh e h)~ X , where J2 P (E/J 
stands for the sum over all particle (hole) lines that 
are cut. The omega line only appears in the sum 
if it is also cut. It enters with a + sign if it is 
directed upwards and with a — sign if it is directed 
downwards. 
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FIG. 8. Basic frequency and finite basis set representation 
PP diagram. 



1 



11(1,2; 3, 4;w) 



+ 




FIG. 9. Basic frequency and real space representation PP 
diagram. 



(6) There is also an overall sign given by the formula 
(— l) h+l , where h is the number of hole lines and I is 
the number of closed loops, including the horizontal 
dotted event lines but ignoring the w lines. 

Diagrams are shown for the independent particle approx- 
imation in Fig. 111! The first diagram reads, 



II a i !a j (a;) — 

The second diagram reads, 
1 



1 



ui + e, - e„ 



-1 



(C3) 



(C4) 



-u> + e l - e a uj + e a - e l 
These two equations are often condensed in the literature 



as. 



lu + e„ 



(C5) 



Let us now introduce one-electron perturbations in the 
form of M circles. 



n(w) 



FIG. 10. Time-unordered representation PP diagram. 
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FIG. 11. Zero-order PP diagrams. 
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p q p q p 

FIG. 12. Electron repulsion integral diagrams. 



q 



(7) Each M circle in a diagram contributes a factor of 
(p\M xc \q), where p is an incoming arrow and q is an 
outgoing arrow. For example, the term correspond- 
ing to Fig. |3] (d) contains a factor of (a\M xc \c), 
while the term corresponding to Fig.|3](j) contains a 
factor of (k\M xc \i). This is a second type of "event" 
(representing "collision" with the quantity M xc ). 

For example, the term corresponding to Fig. [3] (b) is, 



\M xc \b) 



(u> - e k + e c )(efc - e b ) 



(C6) 



This brings us to the slightly more difficult treatment 
of electron repulsions. 

(6) When electron repulsion integrals are represented 
by dotted lines (Feynman and Goldstone dia- 
grams) , each end of the line corresponds to the la- 
bels corresponding to the same spatial point. The 
dotted line representation may be condensed into 
points (Abrikosov and Hugcnholtz diagrams) as in 
Fig. 1121 A point with two incoming arrows, labeled 
p and q, and two outgoing arrows, labeled r and 
s contributes a factor of (pq\\rs) = {pT\fu\qs) — 
(ps\fH\qr). 

(7) To determine the number of loops and hence the 
overall sign of a diagram in which electron repul- 
sion integrals are expanded as dots, then write each 
dot as a dotted line (it does not matter which one of 
the two in Fig.[H]is chosen) and apply rule (6). The 
order of indices in each integral (p(?||rs) should cor- 
respond to the expanded diagrams. (When Gold- 
stone diagrams are interpreted in this way, we call 
them Brandow diagrams.) 
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(8) An additional factor of 1/2 must be added for each Additional information about Hugenholtz and other dia- 
pair of equivalent lines. These are directed lines grams may be found, for example, in Ref. [89j . 
whose interchange, in the absence of further label- 
ing, leaves the Hugenholtz diagram unchanged. 

For example, the term corresponding to Fig. [3] (f) is, 

(ka\\ic) 



(-W + e k - e c )(-CJ + e,i - e a ) 
(ak\\ic) 



(-w + e k - e c )(-w + e, - e a ) 
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